Airfoil Shape and Thickness

Introduction A transient pulse technique is used to The calculationof aeroelastic response in " obtain harmonic forces from a time-marching the transonic regime continues to be an active solution of the complete unsteady transonic area of research because of the significanceof small perturbation potential equation. The flutter in this speed range. Finite-difference unsteady pressures and forces acting on a model codes now allow the analyst to include the nonof the NACA 64A010 conventionalairfoil and the linear aerodynamiceffectsassociatedwith shock MBB A-3 supercritical airfoil over a range of wave formation and motion. These codes have Mach numbers are examined in detail. Flutter been most completely developed for the twocalculationsat constant angle of attack show a dimensional, small-perturbation potential flow similar flutter behavior for both airfoils, equation. The LTRAN2 codeI and its extensions except for a boundary shift in Mach number have been widely used for this purpose. associated with a corresponding Mach number Edwards, et al.2 summarized applications of shift in the unsteady aerodynamic forces, the finite-_-ddif-f-6rence codes to the flutter proDifferences in the static aeroelastic twist blem. In the present paper, a version of the behavior for the two airfoils are significant. LTRAN2 code which incorporateshigher frequency effects3 is used. Nomenclature The aerodynamicforces generatedby a timemarching, finite-differencecode can be incora pitch axis location, referencedto porated into a flutteranalysis in severalways. midchord, in semichords Classical solutions employ the forces for harb semichordlength monic motion in each vibrationmode in an eigenCP center-of-pressurelocation value analysis to determine the flutter point. Cp pressure coefficient The finite-differencecodes may be used to genc chord length erate forces for such an application by forcing c_ lift coefficient the airfoil to oscillate in the prescribed mode c_h lift coefficientdue to plunge at the desired frequency and marching the soluc_ lift coefficientdue to pitch tion in time until the transients have decayed C,n moment coefficientabout c/4 and a periodic solution is achieved. This cmh moment coefficientdue to plunge method requires that the code be run several cm_ moment coefficientdue to pitch times to obtain a single flutter point, a g structuraldamping coefficient second approach is to couple the structuraland h plunge displacement,in semichords aerodynamic equations and time-march the comhl plunge amplitude, in semichords plete system from an initial disturbance. The k reduced frequency,bw/V stability is assessed from the growth or deca_ M free streamMach number of the resulting transient response. Again, m airfoil _ss per unit span several runs of the code are required to deterr_ radius of gyration, referencedto mine the neutral stability point. A third pitch axis, in semichords approach is to compute the harmonic forces from s Laplace transform variable the Fourier transform of the indicial response. t time That is, a step change in the mode is used as v free stream velocity the initialcondition,and the solution is timex streamwisecoordinate marched until the transienthas decayed. xs shock location x_ center-of-gravitylocation,referA transient pulse technique,3 which is a enced to pitch axis, in semichords variation of the indicial response method, is z coordinate normal to stream used in this paper to obtain the aerodynamic ao mean angle of attack, deg forces. The Fourier transform of the resulting al amplitudeof pitch oscillation,deg force response is divided by the transform of , _e static elastic twist angle, deg the airfoilmotion to obtain the harmonic transACp pressure coefficient,normalizedby fer function. A Pad_ fit of these k-plane modal amplitude forces is made to obtain an s-plane representamass ratio, m/_pb2 tion. A root locus solution is then used to p air density determine the stability boundary. This proceT nondimensionaltime, Vt/b dure has the advantage that the expensive aero, @ phase angle, deg dynamic code need be run only once for each Mach oscillationfrequency number, mean angle-of-attack, and vibration _h uncoupled plungemode frequency mode. Subsequent fluttercalculationsfor vari_ uncoupled pitch mode frequency ation_ in structural parameters are relatively inexpensive. The applicabilty of the linear transform techniques to this nonlinear problem described above is verified by comparing the resultswith forced harmonic and with transient aeroelasticsolutions. Edwards, et al.2 have shown a striking Steady Aerodynamics difference in the flutter behavior of a conventional and a supercriticalairfoil in which the The steady state pressure distributionsfor latter shows a curl back in the flutter speed the two airfoils are shown in Fig. 1. The vs. Mach number when static aeroelatic twist is steady angle of attack for each airfoil was included. In this paper, the flutter behavior chosen such that the shock location and lift of a model of the NACA 64A010 conventionalairwere about equal at the same Mach number. Howfoil and the MBB-A3 supercritical airfoil is ever, the shock strength is m_ch less on the examined in detail for one mean angle of attack, aft-loaded supercritical airfoil. Figure 2 The associatedunsteady pressures and forcesare gives the static aerodynamicparameters as funcpresented. Similiarities and differences bettions of Mach number. As noted above, the lift ween the aerodynamics of the two airfoils are coefficients and shock locations are approxiinterpretedto explain the flutter behavior, mately equal, although the lift increases more rapidlywith Mach number on the 64AMES airfoil. As expected, the pitching moments (about the AerodynamicMethod quarter chord) and center-of-pressurelocations are quite different due to the aft loading on All calculations were made with the the supercriticalMBB-A3 airfoil. XTRAN2L3 time-marching finite-difference code which solves the complete unsteady transonic small perturbation potential equation. Th_s Unsteady Aerodynamics code is an enhanced version of the LTRAN2-NLR_ code that includes all of the appropriate timePressure distributions and generalized dependent terms in the differentialequation and aerodynamic forces were computed for harmonic boundary conditions. The alternating-direction oscillationsin two modes, pitch about the quarimplicit algorithm of Rizzetta and Chin5 is ter chord and vertical translation (plunge). used in a three-time-levelscheme to treat the All calculations were made for constant mean second order time derivativein the differential angle of attack. For harmonic oscillation of equation. The low frequency far-field radiation amplitude _1 the pitch motion is given by boundary conditions of Engquist and Majda6 have been extended to the full frequency equa_(T) = _o + _1 sin kT tion and incorporatedinto XTRAN2L. where k = bw/V is the reducedfrequency based on The XTRAN2L code uses a default 80x61 comsemichord and T = Vt/b is the nondimensional putational x-z grid with 51 points on the airtime. The plunge displacementis given by foil. The points are uniformly spaced over the airfoil with the exception of an additional h(T) = hI sin k_ point near the leading edge. The grid extends 20 chord lengths upstream and downstream from where hI is nondimensionalplunge amplitude. the airfoil and 25 chord lengths above and below the airfoil. Reference 3 points out the necesFlutter calculationsgenerally require the sity of stretching the grid smoothly away from determination of generalizedaerodynamic forces the airfoil and of maintaining adequate grid for a range of frequencies for each structural resolution in the far field to avoid erroneous vibration nw)de. Harmonic airloads may be calinternal grid reflection of the outgoing waves, culated with the XTRAN2L code by specifying the Details of the XTRAN2L code and the demonstramode of motion and frequency of oscillation. tion of its accuracy are given in reference 3. Typically four cycles of oscillationare sufficient for the unwanted transients to die out. Fourier analysis of the last cycle of osci_laResults tion gives the Fourier harmonic components of the response. Each such calculationmay require Calculationswere made for two airfoils at as many as 1440 time steps. Because the flutter Fiachnunlbes froi_0.75 to 0.80 in incrementsof analyst is interestedin the forces for a range 0.01. The _irfoils were the _nodelof the NACA of frequencies,the use of an indicial method is 64A010 '_s_ed at NASA Ames Research Center7 attractive. Here the Fourier transform is used (called 64AMES herein) and the supercritical to provide the response at all frequencies of MBB-A3 of German design. These airfoilssare two interest from a single transient response. A of the AGARD standard configurations. The variation of this approach is used herein. 64AMES airfoil is about 10.6 percent thick and has a small amount of camber. The MBB-A3 airPulse Technique. The indicial response is the foil is an 8.9 percent thick, aft-loaded superresponse to a step function change in angle of critical airfoil. The ordinates of both airattack. For the equation solved by the XTRAN2L foils are given in reference 8. The airfoil code with the second order time derivative slopes required by the aerodynamic code were included, the resulting response includes nonobtained from cubic spline fits to the airfoil physical transients due to the approximation of abscissas and ordinates, using an approximation the infinite initial derivative using the to airfoil arc length as the independent finite time steps taken. For this reason variable. All calculationswere made at mean smoothly varying exponentially shaped pulsej angles of attack of _o = 10 for the 64AMES is used. The input pulse in angle of attack is airfoil and _o = -0.50 for the MBB-A3 airgiven by foil. _(T) = _0 + _Ie-(T-17.5_T)2/4 where L_r is the nondimensionaltime step (_I = obtained from forced harmonic oscillationof the 0.50 and A_ = 5_/32 were used), airfoil at k = 0.15. The pressure pulse seen near 55% chord is a result of the shock wave An example of the input pulse and resulting motion. The greatest nonharmonic response lift and moment response for the 64AMES airfoil occurred at the aft end of the shock pulse and at M = 0.8 is shown in Fig. 3. The forces have was as large as the fundamental component for not returned to their initialvalues for the 128 _1 = 20. Away from the shock, the response time steps shown. It is essential that the calwas nearly all at the fundamental,frequency. culation be carried out to sufficient time for the transient to decay completely for the A phase jump of about 1200 occurs at the following reason. The harmonic response is aft end of the shock pulse. Away from the obtained from the fast Fourier transforms (FFT) shock, the phases and normalized magnitudes are of the input and the response. If the final essentially equal. The width of the pulse value of the time history has not returned to increases while the height decreases with the initial value, the low frequency results increasing _1 in such a way that the overall will be seriously in error. In practice 2048 forces remain nearly equal. That is, the force steps were used, although 1024 steps were somecoefficientsare independentof the amplitude of times sufficient. The FFT technique provides oscillation. The maximum differencesin magniforces up to the Nyquist frequency of tude and phase for the cases shown are about k = :r/AT= 6.4. The accuracy of the method for 0.5% and 10 for lift and 6.4% and 70 for frequencies as high as k = 2 is demonstratedin moment. reference3. Harmonic Airloads. Calculations for Steady results may be obtained from the forced harmonic motion were made for the plunge XTRAN2L code using either the steady, successive and pitch modes at each Mach number for both line over-relaxation solver, or the unsteady airfoils at k = 0.15. The plunge amplitude was solver with fixed airfoil geometry and timehl = 0.02, and the pitch amplitude was marching to a steady state. These options lead _1 = 0.50. The magnitude and phase of the to different steady state results as indicated lifting pressure coefficients (normalized by in Fig. 4 for the 64AMES airfoil at M = 0.8. modal amplitude) for these cases are shown in Here the steady solver was run to convergence Fig. 7 for the pitch mode and Fig. 8 for the with a lift coefficient c_ = 0.3461 obtained, plunge mode. The chief difference between the The unsteady solverwas then used (with the airtwo airfoils is in the strength of the shock foil held fixed) which resulted in the transient pulses; the location is about the same. This shown in Fig. 4. As was the case for the pulse result might be anticipatedfrom the corresponddescribed in the preceding paragraph, it was ing steady pressures shown in Fig. 1. At the necessary to calculate several thousand time same Mach number, the phase jump at the shock is steps for this transient to decay. The transomewhat less for the MBB-A3 airfoil than for sient shown was for the first 512 time steps, the 64AMES airfoil for the same reason. The resulting lift coefficient at convergence was c_ = 0.3354. It is, therefore, necessary The unsteady forces that result from these to run the unsteady solver with fixed airfoil pressure distributionsat k = 0.15 are shown in before applying the pulse technique described Figs. 9 and 10. Linear theory resultsfrom kerabove in order to insure that the initial and nel function calculationsfor a flat plate airfinal conditions for the pulse are the same. foil are also shown. The lift coefficients (Fig. 9) for both airfoils show only mild Harmonic Forces. The harmonic forceswere changes with Mach number and agree quite well. obtained using the pulse techniqueby performing The moment coefficients do show significant an FF-Tof the pulse input (Fig. 3(a)) and of the changes with Mach number. The most interesti_g force output (Fig. 3(b)) and dividing the latter observationfrom these two figures is the agreeby the former. For a case with 2048 time steps ment between the two sets of calculationswhen this leads to a frequency resolution of the results for the MBB-A3 airfoil are shifted Ak = 1/160. A typical result is shown in Fig. to a 0.01 lower F1ach number. That is, the 5. The real and imaginary components,of the forces for the two airfoils are nearly equal four force coefficientsare presented for t_'oacross the _ach number range if the square symcases the 64AMES airfoi! at M = 0.78 and the bols are shifted 0.01 Mach number to the left. MBB-A3 airfoil at M = 0.79. The most striking This conclusionwas anticipatedin Fig. 5, which result is the agreement in the forces for these shows that this result holds over a wide fretwo different airfoils at the different Mach quency range. numbers. This agreement in the forces with a 0.01 shift in M holds true across the Mach numIn order to verify the accuracy of the bet range for the cases treated in this paper, pulse technique,the forces shown in Figs. 9 and AS may be anticipated, this comparisonwill be 10 were comparedwith those obtained from pulses reflected in the flutter results to be shown in h and _. The maximum differences in magnilater. The jagged nature of the curves at low tude occurred for the pitch case and were frequency is not a serious problem in flutter 0.048 for c_ and 0.018 for The analysis for which k is usually greater than maximum phase _fference was 30. Cma" 0.1. OscillationAmplitude Effects. The liftFlutter in9 pressure distributionon the 64AMES airfoil at M = 0.78 oscillatingin pitch at four differFlutter boundarieswere calculated for each ent amplitudes is shown in Fig. 6. The first airfoil at each Mach number and one mean angle harmonic only is shown. These results were of attack (_o = 1o for the 64AMES airfoil coincident. This result is to be expected in and ao = -0.5° for the MBB-A3 airfoil). The light of the nature of the aerodynamic forces structural parameters were those of Isogai's shown earlier (Figs. 5, 9, and 10). The simiCase A9 which were used in reference 2. The larity in the flutter boundarieswas surprising elastic-axis location was at a = -2, the center because of the differences in flutter behavior of gravity at x_ = 1.8, and the radius of for the two airfoils reported in reference 2. gyration squared was r_2 = 3.48. The mass Additional calculations made to explain these ratio was u = 60 and the modal frequencieswere differencesare described in the.followingpara_h = _ : 100 rad/sec, graphs. Flutter Solution Method. The XTRAN2L code Figure 13 shows the present flutter boundmay be used to obtain aeroelastic solutions by aries as lines with the results from reference 2 time-marchingthe coupledaerodynamicand strucand some additional calculations shown as symtural equations.2 This method requires severbols. The flagged-symbol points, taken from al runs of the code to determine the neutral reference 2, were calculated using the LTRAN2stability (flutter) point. Alternatively, the NLR code4 and differ from the present XTRAN2L harmonic forces obtained from the pulse technicode results in: (I) equationlevel (the second que described herein may be used in a conventime derivative term in the governing equation tional V-g flutter analysis. The harmonic was omitted), (2) finite-difference grid forces (k-plane) may also be representedin the (discussed in reference 3), and (3) airfoil Laplace variable (s-plane), and root locus slope calculation procedure. In addition, the methods may be used to determine the flutter NACA 64A010 was used rather than the 64AMES airbehavior. In this paper a Pad_ fit10 was used foil used herein. In an attempt to sort out to obtain the s-plane representation. It was these different effects, several new calculanecessary to use care in applying this method tions were made. because of inaccuraciesin the k-plane forces at very low frequencies (k<O.1). There is also The unflaggedcircles shown in Fig. 13 were some question as to the suitability of the computed for the NACA 64A010 airfoil using the particular form of fit for transonic airloads. XTRAN2L code. The difference between these In spite of these reservations, the method was points and the solid curve for the 64AMES airused successfullyby insuring that the Pad_ fit foil is a thickness effect and is consistent was accurate near the expected flutter frewith a transonicsimilarityshift of about 0.007 quency, in M corresponding to a ratio in thickness of 1.06 between the airfoils. The difference betA sample ro()tlocus flutter solution for ween the flaggedand unflagged circles (both for the 64AMES airfoil is shown in Fig. 11. The the 64A010 airfoil) is due to the three code different curves were obtained by varying the differencesenumeratedabove. speed. The symbols represent the time-marching aeroelastic solutions. The tic marks on the The unflagged square symbols for the MBB-A3 curves occur dt speeds which correspondto those airfoil, shown in Fig. 13, were obtained from for the symbols. Excellent agreement is shown, the XTRAN2L code with the second time derivative which verifies the chain of linearity assumpterm set to zero. The rather large difference tions used in obtaining the curves, i.e. pulse between these symbols and the dashed curve for transient to harmonic forces (via FFT) to root the MBB-A3 airfoil is due solely to the locus (via Pad_ fit). difference in equation level (item (1) above). The flagged and unflagged squares were computed Flutter Results. The flutter boundaries using the same equation level and differ only in for the two airfoils are shown i:;FiL_.12 as items (2) and (3) above. This differenceis not flutter speed index V/b_V_ an(! flutter large, but the results from reference 2 reduced frequency as functions of Mach number. (flagged) do show a steeper transonic dip in The Mach number range shown covers the region of flutter speed. In particular, at M = 0.8 the the transonic dip but does not extend past the flutter speed from reference 2 (flagged square) point of minimum flutter speed. The potential has dropped below the speed computed from the flow calculationsbecome unreliable beyond about same equation (unflagged square) herein. This M = 0.8 for the 64AMES airfoil. Difficulties difference in trend emphasizes the sensitivity with convergence of the aerodynamic code occur of transonic calculations to the finitewhich appear to correlatewith known problems of differencegrid and airfoil slopes used. nonuniqueness with potential flow theory.11 An example of this problem was observed in a Static Aeroelastic Twist. The flutter calculation for the 64AMES airfoil at M = 0.84 resultsshown in Figs. 12 and 13 were calculated and ao = OO for which a decidedly nonsymholding the mean angle of attack of the section metric pressure distributionwas obtained; the fixed. From a strip theory point of view, this upper and lower surface shocks were located at would require a different root trim angle for x/c = 0.82 and 0.69, respectively. For this each point because of the variations in flutter airfoil at M = 0.81 and _o = lO (the angle speed and in static pitching moment coefficient of attack for the flutter calculations) the at each point. The elastic twist is computed upper surface shock was located at x/c = 0.82 from after 2048 time steps and the solution was not 2 V )2Cm(%,M ) converged, a_e = _o_r = _ ( irr The flutter boundariesfor the two airfoils (Fig. 12) are remarkablysimilar. If the boundwhere _o is the section angle of attack and ary for the MBB-A3 airfoil is shifted to the _r is the unknown root angle required to left by 0.01 Mach number, the curves are nearly balance the static aeroelastic moment. The pitching moment coefficient Cm(_o,M) depends 3. Flutter boundaries were calculated for a on _o (_o = I° and -0.5o for the 64AMES conventionaland a supercrltical airfoil. The and MBB-A3 airfoils, respectively) and on Mach mean angles of attack were chosen such that the number (as shown in Fig. 2 for cm referenced steady lift and shock locations were the same. to c/4) and is referenced to the twist axis The flutter boundarieswere veryslmilar except here. for a shift of about 0.01 in Mach number. , Reference 2 showed a marked difference in 4. The unsteady forces for the two airfoils behavior for the flutterof the conventionaland were compared at a reduced frequency (k=0.15) in supercriticalairfoilswith Mach numberwhen the the flutter range. The lift coefficientswere root trim angle was held fixed. (See Figs. 15 quite similarbut the moment coefficientsshowed , and 16 of reference 2 in which the _o corresome differences. Again, a shift of 0.01 in sponds to _r herein). The determination of Mach number brought the forces into reasonable those results required calculation for several agreement. values of _o and interpolation of the resulting root angles to obtain the desired ar. In 5. The unsteady pressure distributions for the this paper, all calculationswere made for contwo airfoils were also compared at a reduced stant airfoil section angle _o. Using the frequencyof 0.15. Although the shock strengths root locus procedure, it was easy to vary the were different, the locationsand widths of the section elastic axis location and the result of shock pulses were very similar. this variation on the static elastic twist is describedbelow. 6. Unsteady pressure distributionsfor a range of pitch oscillation amplitudes showed that In Fig. 14 the elastic twist at flutter is although the width of the shock pulse increased shown for each airfoil for three elastic axis with increased amplitude, the unsteady coeffilocations. The value a = -2 was used in all of cients were nearly constant. the results shown earlier and the twist angles for these curves correspond to the flutter 7. Calculationsof the static aeroelastictwist boundaries of Fig. 12. Comparing the two airat flutter were made for several pitch axis lofoils for a = -2, the MBB-A3 airfoil has a larcations. The supercritical airfoil showed a ger twist and shows a greaterMach number effect greater variation of twist with Mach number than than does the 64AMES airfoil. The aeroelastic did the conventionalairfoil. windup with M of the MBB-A3 airfoil (2.50) is about one degree greater than that of the 64AMES References airfoil (1.5o). This greater twist is a result of the larger pitching nDment coefficient IBallhaus, W. F. and Goorjian, P. M.: for the supercriticalairfoil (Fig. 2). If calImplicit Finite-Difference Computations of culations were made for constant _r one would Unsteady Transonic Flow About Airfoils. AIAA expect the flutter similiarities shown in Fig. Journal, vol. 15, no. 12, Dec. 1977, pp. I"T_I_-_12 to disappear. The shapes of the flutter boundaries with Mach number would be quite different because the airfoil angle of attack _o 2Edwards, J. W., Bennett, R. M., Whitlow, would change i_re for the MBB-A3 than for the W., Jr., and Seidel, D. A.: Time-MarchingTran64AMES airfoil, sonic Flutter Solutions Including Angle-ofAttack Effects. AIAA paper No. 82-0685, presentThe value a = -2 places the pitch mode axis ed at the AIAA/ASME/ASCE/AHS 23rd Structures, one semichord length ahead of the leadingedge. Structural Dynamics and Materials Conference, As this axis is moved closer to the lea,lingedge New Orleans, Louisiana,May 1982. (a = -1.8 and -1.6), Fig. 14 shows the resulting decrease in elastic twist and a less severe Mach 3Seidel, D. A., Bennett, R. M., and number trend. At a = -1.6, the two airfoils Whitlow, W., Jr.: An Exploratory Study of show nearly the sa_ Mach number effect. For Finite-DifferenceGrids for Transonic Unsteady this variation the airfoil mass propertieswere Aerodynamics. AIAA Paper No. 83-0503, presented held fixed, that is a + x_ = -0.2 in all at the AIAA 21st Aerospace Sciences Meeting and cases. Technical Display, Reno, Nevada,January 1983. 4Houwink, R. and van der Vooren, J.: Conclusions Improved Version of LTRAN2 for Unsteady Transonic Flow Computations. AIAA Journal, vol. 18, 1. A transient pulse technique has been shown no. 8, Aug. 1980, pp. 1008-1010. o to be an accurate and efficient tool for determining aerodynamic forces for use in transonic 5Rizzetta, D. P. and Chin, W. C.: Effect flutter analysis. The necessity of carefully of Frequency in Unsteady Transonic Flow. AIAA monitoring the convergenceof the time-marching Journal, vol. 17, no. 7, July 1979, pp. 779-'TB_I_ -. calculation is emphasized. 6Engquist, B. and Majda, A.: Numerical 2. An s-plane Pad6 representationfor the harRadiation Boundary Conditions for Unsteady Tranmonic forces was used to provide a convenient sonic Flow. Journal of Computational Physics, and efficient means for performing parameter vol. 40, 1981, pp. 91-103. variations in flutter analysis. The technique used was sensitive to small inaccuraciesin the 7Davis, Sanford S. and Malcolm, Gerald harmonic forces at low frequencies. N.: ExperimentalUnsteady Aerodynamicsof Conventional and SupercriticalAirfoils. NASA TM 81221, August 1980. 8Bland, S. R.: AGARD Two-Dimensional 1OEdwards, O. W.: Applications of Aeroelasttc Configurations. AGARD Advisory Laplace Transform Methods to Airfoil Motion and Report No. 156, August 1979. Stability Calculations. AIM Paper No. 79-07721 AIAA Structures, Structural Dynamics, and 9Isogai, K.: Numerical Study of TranMaterials Conference, St. Louis, Ntssourf, May sonic Flutter of a Two-Dimensional Airfoil. 1979. Nattonal Aerospace Laboratory, Japan, TR-617T, llStetnhoff, John and Jameson, Antony: July 1980. Multiple Solutions of the Transonic Potential Flow Equation. AIM paper No. 81-1019, AIAA Computational Fluid Dynamics Conference, Palo UPPER Alto, California, June 1981 S RFACE M= 0.80 l.O .79 /" UPPER M = 0.80 .18 /SURFACE jr'j/-. 79 37 .78 -.5 cm 0 .2 .4 .6 .8 l.O 0 .2 .4 .6 .8 I.O